This notebook pulls in data on Nipah RBP entry into CHO-EFNB2 and CHO-EFNB3 cells, filters data, calculates stats, and makes figures¶

In [1]:
# this cell is tagged as parameters for `papermill` parameterization
e2_distances_file = None
func_scores_E2_file = None
func_scores_E3_file = None

filtered_E2_data = None
filtered_E3_data = None
contact_type_plot = None
E2_E3_entry_corr_plot = None
E2_E3_entry_all_muts_plot = None
combined_E2_E3_correlation_plots = None
entry_region_boxplot_plot = None

nipah_config = None
altair_config = None

entropy_file = None
surface = None
In [2]:
# Parameters
func_scores_E2_file = "results/func_effects/averages/CHO_EFNB2_low_func_effects.csv"
func_scores_E3_file = "results/func_effects/averages/CHO_EFNB3_low_func_effects.csv"
e2_distances_file = "results/distances/2vsm_distances.csv"
contact_type_plot = "results/images/contact_type_plot.html"
filtered_E2_data = "results/filtered_data/E2_entry_filtered.csv"
filtered_E3_data = "results/filtered_data/E3_entry_filtered.csv"
E2_E3_entry_corr_plot = "results/images/E2_E3_entry_corr_plot.html"
E2_E3_entry_all_muts_plot = "results/images/E2_E3_entry_all_muts_plot.html"
combined_E2_E3_correlation_plots = (
    "results/images/combined_E2_E3_correlation_plots.html"
)
nipah_config = "nipah_config.yaml"
altair_config = "data/custom_analyses_data/theme.py"
entropy_file = "results/entropy/entropy.csv"
entry_region_boxplot_plot = "results/images/entry_region_boxplot_plot.html"
surface = "data/custom_analyses_data/surface_exposure.csv"
In [3]:
import math
import os
import re
import altair as alt

import numpy as np

import pandas as pd

import scipy.stats

import Bio.SeqIO

import yaml

from Bio import AlignIO
from Bio import PDB
from Bio.Align import PairwiseAligner
In [4]:
# allow more rows for Altair
_ = alt.data_transformers.disable_max_rows()

if os.getcwd() == '/fh/fast/bloom_j/computational_notebooks/blarsen/2023/Nipah_Malaysia_RBP_DMS/':
    pass
    print("Already in correct directory")
else:
    os.chdir("/fh/fast/bloom_j/computational_notebooks/blarsen/2023/Nipah_Malaysia_RBP_DMS/")
    print("Setup in correct directory")
Setup in correct directory
In [5]:
if surface is None:
    e2_distances_file = "results/distances/2vsm_distances.csv"
    func_scores_E2_file = "results/func_effects/averages/CHO_EFNB2_low_func_effects.csv"
    func_scores_E3_file = "results/func_effects/averages/CHO_EFNB3_low_func_effects.csv"
    
    filtered_E2_data = "results/filtered_data/E2_entry_filtered.csv"
    filtered_E3_data = "results/filtered_data/E3_entry_filtered.csv"
    
    nipah_config = "nipah_config.yaml"
    altair_config = "data/custom_analyses_data/theme.py"
    
    entropy_file = "results/entropy/entropy.csv"
    surface = "data/custom_analyses_data/surface_exposure.csv"

Read in custom altair theme and Import YAML file with parameters¶

In [6]:
if altair_config:
    with open(altair_config, 'r') as file:
        exec(file.read())

with open(nipah_config) as f:
    config = yaml.safe_load(f)

Filter and merge EFNB2 and EFNB3 entry¶

In [7]:
#Import median entry scores from different selections
func_scores_E2 = pd.read_csv(func_scores_E2_file).dropna().round(3) # this removes wildtype observations from df
func_scores_E3 = pd.read_csv(func_scores_E3_file).dropna().round(3) # this removes wildtype observations from df
In [8]:
def num_selections_and_filter(df, name):
    def calculate_filtering_mutants(df, name):
        print(f"\nThe dataset is: {name}")
        
        # How many selections were performed
        max_sels = df["n_selections"].max()
        num_sels_cutoff = (max_sels / 2) + 1

        # Find total number of possible mutants
        total_mut = (602 - 71) * 19
        
        # Filter data
        filter_test = df[
            (df["site"] != 603) & (df["mutant"] != "-") & (df["mutant"] != "*")
        ]

        num_variants_pre_filter = filter_test.shape[0]
        print(
            f"After filtering stop and gaps, there are {num_variants_pre_filter} mutants which is {(num_variants_pre_filter/total_mut) * 100:.1f}%"
        )

        filter_test_times_seen = filter_test[
            filter_test["times_seen"] >= config["func_times_seen_cutoff"]
        ]
        num_variants_times_seen = filter_test_times_seen.shape[0]
        print(
            f"After filtering for {config['func_times_seen_cutoff']} times seen, there are {num_variants_times_seen}, which is {(num_variants_times_seen/total_mut) * 100:.1f}%"
        )

        filter_test_effect_std = filter_test_times_seen[
            filter_test_times_seen["effect_std"] <= config["func_std_cutoff"]
        ]
        num_variants_std = filter_test_effect_std.shape[0]
        print(
            f"After filtering for {config['func_std_cutoff']} std cutoff, there are {num_variants_std}, which is {(num_variants_std/total_mut) * 100 :.1f}%"
        )

        filter_test_n_selections = filter_test_effect_std[
            filter_test_effect_std["n_selections"] >= num_sels_cutoff
        ]
        num_variants_n_selections = filter_test_n_selections.shape[0]
        print(
            f"After filtering for mutants in in all selections, there are {num_variants_n_selections}, which is {(num_variants_n_selections/total_mut) * 100 :.1f}%"
        )

    def apply_filters(df, name):
        # Now do the main filtering
        max_sels = df["n_selections"].max()
        num_sels_cutoff = (max_sels / 2) + 1
        print(
            f"The number of selections a mutant must be observed in is: {num_sels_cutoff}"
        )
        # The main filtering. Filters site 603 (is a stop codon/end of gene and we don't want those mutants). Also filter out stop mutants and apply filtering from config file
        filtered_df = df[
            (df["site"] != 603)
            & (df["mutant"] != "-")
            & (df["mutant"] != "*")
            & (df["times_seen"] >= config["func_times_seen_cutoff"])
            & (df["effect_std"] <= config["func_std_cutoff"])
            & (df["n_selections"] >= num_sels_cutoff)
        ]
        return filtered_df

    # Filtering stats
    calculate_filtering_mutants(df, name)  # call definition above

    # Filter
    filtered_df = apply_filters(df, name)

    # Now write filtered_data to results/
    if name == "E2":
        filtered_df.to_csv(filtered_E2_data, index=False)
    if name == "E3":
        filtered_df.to_csv(filtered_E3_data, index=False)

    # return filtered dataframe
    return filtered_df

# Call the filtering functions
func_scores_E2 = num_selections_and_filter(func_scores_E2, "E2")
func_scores_E3 = num_selections_and_filter(func_scores_E3, "E3")

# make a merged dataframe of ephrin-b2 and ephrin-b3 entry data
def merge_data(df1,df2):
    merged_df = pd.merge(
        df1,
        df2,
        on=["site", "mutant", "wildtype"],
        how="outer",
        suffixes=["_E2", "_E3"],
    )
    df1["cell_type"] = "CHO-EFNB2"
    df2["cell_type"] = "CHO-EFNB3"
    concat_df = pd.concat([df1,df2])

    return merged_df,concat_df

merged_df,concat_df = merge_data(func_scores_E2,func_scores_E3)

# Show some stats of filtered merged data
stats = merged_df.describe().round(2)
display(stats)
The dataset is: E2
After filtering stop and gaps, there are 10040 mutants which is 99.5%
After filtering for 2 times seen, there are 9789, which is 97.0%
After filtering for 1 std cutoff, there are 9729, which is 96.4%
After filtering for mutants in in all selections, there are 9725, which is 96.4%
The number of selections a mutant must be observed in is: 5.0

The dataset is: E3
After filtering stop and gaps, there are 10072 mutants which is 99.8%
After filtering for 2 times seen, there are 9852, which is 97.7%
After filtering for 1 std cutoff, there are 9714, which is 96.3%
After filtering for mutants in in all selections, there are 9586, which is 95.0%
The number of selections a mutant must be observed in is: 4.5
site effect_E2 effect_std_E2 times_seen_E2 n_selections_E2 effect_E3 effect_std_E3 times_seen_E3 n_selections_E3
count 9885.00 9725.00 9725.00 9725.00 9725.00 9586.00 9586.00 9586.00 9586.00
mean 337.22 -0.92 0.32 7.20 7.99 -0.95 0.34 6.14 6.98
std 153.23 1.40 0.24 4.08 0.12 1.46 0.25 3.27 0.14
min 71.00 -3.52 0.00 2.00 5.00 -3.61 0.00 2.00 5.00
25% 205.00 -2.21 0.15 4.62 8.00 -2.50 0.17 4.14 7.00
50% 338.00 -0.27 0.30 6.38 8.00 -0.22 0.33 5.43 7.00
75% 470.00 0.22 0.47 8.50 8.00 0.22 0.49 7.14 7.00
max 602.00 0.64 1.00 64.38 8.00 0.66 1.00 49.14 7.00

Stats¶

In [9]:
def calculate_stats(df, name):
    print(f"For {name}:")
    total_mut = (602 - 71) * 19
    print(f'There are {total_mut} amino acid mutations possible')
    muts_present = df["effect"].shape[0]
    fraction_muts = muts_present / total_mut
    print(
        f"fraction muts present in {name} is {fraction_muts:.2f} {muts_present}/{total_mut}"
    )

    deleterious_muts = df[df["effect"] <= -0.25].shape[0]
    neutral_muts = df[(df["effect"] > -0.25) & (df["effect"] < 0.25)].shape[0]
    positive_muts = df[df["effect"] > 0.25].shape[0]

    frac_bad_muts = deleterious_muts / muts_present
    frac_neutral_muts = neutral_muts / muts_present
    frac_pos_muts = positive_muts / muts_present
    print(
        f"The number of deleterious mutants for {name} is {frac_bad_muts:.2f} {deleterious_muts}/{muts_present}"
    )
    print(
        f"The number of neutral mutants for {name} is {frac_neutral_muts:.2f} {neutral_muts}/{muts_present}"
    )
    print(
        f"The number of positive mutants for {name} is {frac_pos_muts:.2f} {positive_muts}/{muts_present}\n"
    )

calculate_stats(func_scores_E2, "E2")
calculate_stats(func_scores_E3, "E3")
For E2:
There are 10089 amino acid mutations possible
fraction muts present in E2 is 0.96 9725/10089
The number of deleterious mutants for E2 is 0.50 4907/9725
The number of neutral mutants for E2 is 0.26 2571/9725
The number of positive mutants for E2 is 0.23 2239/9725

For E3:
There are 10089 amino acid mutations possible
fraction muts present in E3 is 0.95 9586/10089
The number of deleterious mutants for E3 is 0.49 4722/9586
The number of neutral mutants for E3 is 0.28 2680/9586
The number of positive mutants for E3 is 0.23 2178/9586

How many sites and which sites only have negative entry scores for mutations?¶

In [10]:
def overall_stats_all_neg(df,effect):
    filtered_df = df.groupby('site').filter(lambda group: (group[effect] < 0).all())
    unique = filtered_df['site'].unique()
    print(list(unique))
    total_sites = df['site'].unique().shape[0]
    subset = filtered_df['site'].unique().shape[0]
       
    fraction = subset/total_sites
    percent = fraction * 100
    print(f'The total number of sites are: {total_sites}')
    print(f' The number of sites where all mutants are negative for {effect}: {subset}')
    print(f' The percent of sites where all mutants are negative for {effect}: {percent:.0f}')
    return unique

intolerant_sites_E2 = list(overall_stats_all_neg(func_scores_E2,'effect'))
intolerant_sites_E3 = list(overall_stats_all_neg(func_scores_E3,'effect'))
[80, 106, 107, 108, 111, 112, 113, 120, 121, 125, 126, 127, 130, 138, 146, 151, 157, 158, 159, 162, 163, 165, 167, 172, 189, 203, 205, 206, 207, 208, 216, 229, 240, 246, 251, 253, 254, 257, 258, 259, 260, 262, 263, 264, 266, 267, 303, 323, 331, 347, 355, 382, 387, 395, 412, 460, 467, 487, 489, 493, 499, 500, 503, 506, 537, 563, 565, 574, 594, 598]
The total number of sites are: 532
 The number of sites where all mutants are negative for effect: 70
 The percent of sites where all mutants are negative for effect: 13
[95, 100, 106, 107, 108, 111, 112, 113, 121, 125, 126, 138, 146, 158, 162, 163, 201, 203, 206, 207, 216, 229, 240, 243, 248, 251, 253, 257, 258, 259, 260, 263, 266, 303, 347, 352, 355, 368, 382, 387, 389, 395, 412, 439, 458, 460, 467, 486, 487, 489, 493, 494, 497, 499, 500, 501, 503, 504, 505, 506, 510, 526, 531, 532, 533, 537, 556, 557, 563, 565, 573, 574, 579, 581, 584, 585, 588, 594]
The total number of sites are: 532
 The number of sites where all mutants are negative for effect: 78
 The percent of sites where all mutants are negative for effect: 15
In [11]:
def calculate_top_func_scores(df,effect):
    percentile_95_effect_E2 = df[effect].quantile(0.999)
    cutoff_E2_df_sites = df[df[effect] > percentile_95_effect_E2]
    E2_site_cutoff = cutoff_E2_df_sites['site'].unique()
    print(f'The sites with the highest functional scores for {effect} are: {list(E2_site_cutoff)}')

calculate_top_func_scores(func_scores_E2,'effect')
calculate_top_func_scores(func_scores_E3,'effect')
The sites with the highest functional scores for effect are: [280, 305, 407, 463, 468, 501, 552, 556, 584, 599]
The sites with the highest functional scores for effect are: [183, 315, 337, 358, 378, 393, 418, 419, 554, 597]

Make bubble plots of receptor contact site type (median values per site)¶

In [12]:
def make_bubbleplot_entry_region(df):  # Create a bubble plot using Altair for contact site mutants
    barrel_ranges = {
        "Hydrophobic": config["hydrophobic"],
        "Salt Bridges": config["salt_bridges"],
        "Hydrogen Bonds": config["h_bond_total"],
        "Contact": config["contact_sites"],
        "Overall": list(range(71, 602)),
    }
    custom_order = [
        "Hydrophobic",
        "Salt Bridges",
        "Hydrogen Bonds",
        "Contact",
        "Overall",
    ]
    empty_chart = []
    variant_selector = alt.selection_point(
        on="mouseover", empty=False, fields=["site"], value=1
    )
    for selection in ["CHO-EFNB2", "CHO-EFNB3"]:
        agg_means = []
        tmp_df = df[df["cell_type"] == selection]
        mean_df = tmp_df.groupby("site")[["effect"]].median().reset_index()

        # For each barrel, filter the site_means dataframe to the sites belonging to that barrel and then store the means
        for barrel, sites in barrel_ranges.items():
            subset = mean_df[mean_df["site"].isin(sites)]
            for _, row in subset.iterrows():
                agg_means.append(
                    {"barrel": barrel, "effect": row["effect"], "site": row["site"]}
                )
        agg_means_df = pd.DataFrame(agg_means)
        chart = (
            alt.Chart(agg_means_df, title=f"{selection}")
            .mark_point()
            .encode(
                x=alt.X(
                    "barrel:O",
                    sort=custom_order,
                    title='Contact Type',
                    axis=alt.Axis(labelAngle=-90),
                ),
                y=alt.Y(
                    "effect",
                    title="Median Cell Entry",
                    axis=alt.Axis(grid=True, tickCount=4),
                ),
                xOffset="random:Q",
                tooltip=["barrel", "effect", "site"],
                size=alt.condition(
                    variant_selector, alt.value(100), alt.value(25)
                ),
                color=alt.condition(
                    variant_selector, alt.value("orange"), alt.value("black")
                ),
                opacity=alt.condition(variant_selector, alt.value(1), alt.value(0.3)),
            ).properties(width=config['width'],height=config['height'])
            .transform_calculate(random="sqrt(-1*log(random()))*cos(2*PI*random())")
        )
        empty_chart.append(chart)
    combined_effect_chart = (
        alt.hconcat(*empty_chart)
        .resolve_scale(y="shared", x="shared", color="independent")
        .add_params(variant_selector)
    )
    return combined_effect_chart


tmp_img = make_bubbleplot_entry_region(concat_df)
tmp_img.display()
if entry_region_boxplot_plot is not None:
    tmp_img.save(contact_type_plot)

Make bubble plots depending on amino acid property¶

In [13]:
def make_bubbleplot_wildtype_prop(df):
    barrel_ranges = {
        "hydrophobic": list(["A", "V", "L", "I", "M"]),
        "aromatic": list(["Y", "W", "F"]),
        "positive": list(["K", "R", "H"]),
        "negative": list(["E", "D"]),
        "hydrophilic": list(["S", "T", "N", "Q"]),
        "special": list(["C", "P", "G"]),
    }
    empty_charts = []
    variant_selector = alt.selection_point(
        on="mouseover", empty=False, nearest=True, fields=["site"], value=1
    )
    for selection in ["CHO-EFNB2", "CHO-EFNB3"]:
        if selection == "CHO-EFNB2":
            effect_name = "EFNB2"
        else:
            effect_name = "EFNB3"

        tmp_df = df[df["cell_type"] == selection]

        unique_wildtype_df = tmp_df[["site", "wildtype"]].drop_duplicates()
        mean_df = tmp_df.groupby("site")[["effect"]].median().reset_index()
        mean_df = pd.merge(mean_df, unique_wildtype_df, on="site", how="left")

        agg_means = []

        # For each barrel, filter the site_means dataframe to the sites belonging to that barrel and then store the means
        for barrel, sites in barrel_ranges.items():
            subset = mean_df[mean_df["wildtype"].isin(sites)]
            for _, row in subset.iterrows():
                agg_means.append(
                    {"wildtype_class": barrel, "effect": row["effect"], "site": row["site"], "wildtype": row["wildtype"]}
                )
        agg_means_df = pd.DataFrame(agg_means)

        chart = (
            alt.Chart(agg_means_df, title=f"{selection}")
            .mark_point()
            .encode(
                x=alt.X(
                    "wildtype_class:O",
                    title="Wildtype amino acid property",
                    axis=alt.Axis(labelAngle=-90),
                ),  
                y=alt.Y(
                    "effect",
                    title=f"Median Cell Entry",
                    axis=alt.Axis(grid=True, tickCount=4),
                ),
                xOffset="random:Q",
                tooltip=["wildtype_class", "effect", "site","wildtype"],
                size=alt.condition(variant_selector, alt.value(100), alt.value(25)),
                color=alt.condition(
                    variant_selector, alt.value("orange"), alt.value("black")
                ),
                opacity=alt.condition(variant_selector, alt.value(1), alt.value(0.2)),
            ).properties(width=config['width'],height=config['height'])
            .transform_calculate(random="sqrt(-1*log(random()))*cos(2*PI*random())")
            
        )
        empty_charts.append(chart)
    combined_effect_chart = (
        alt.hconcat(*empty_charts)
        .resolve_scale(y="shared", x="shared", color="independent")
        .add_params(variant_selector)
    )
    return combined_effect_chart


wildtype_aa_bubble_img = make_bubbleplot_wildtype_prop(concat_df)
wildtype_aa_bubble_img.display()

Plot correlations between E2 and E3 entry¶

In [14]:
# Import distance data
e2_distances = pd.read_csv(e2_distances_file)
distance_df = pd.merge(
    merged_df, e2_distances[["site", "distance"]], on="site", how="left"
)


def determine_distance(df):
    # Define the conditions
    conditions = [
        df["distance"] < 4,
        (df["distance"] >= 4) & (df["distance"] <= 8),
        df["distance"] > 8,
    ]

    # Define the associated values for the conditions
    choices = ["contact", "close", "distant"]

    # Apply the conditions and choices to the 'E2_contact' column
    df["contact"] = np.select(conditions, choices, default="distant")
    return df


distance_df = determine_distance(distance_df)


def median_correlation_plot(df, metric):
    aggregation = getattr(df.groupby("site")[["effect_E2", "effect_E3"]], metric)
    means = aggregation().reset_index()
    slope, intercept, r_value, p_value, std_err = scipy.stats.linregress(
        means["effect_E2"], means["effect_E3"]
    )
    display(r_value.round(2))

    means = means.rename(
        columns={"effect_E2": f"{metric}_E2", "effect_E3": f"{metric}_E3"}
    )

    contact_sites = df[["site", "contact", "wildtype"]].drop_duplicates()
    df_mean = pd.merge(means, contact_sites, on="site", how="left")

    chart = (
        alt.Chart(df_mean)
        .mark_point()
        .encode(
            x=alt.X(f"{metric}_E2", title="Summed Cell Entry for EFNB2"),
            y=alt.Y(f"{metric}_E3", title="Summed Cell Entry for EFNB3"),
            tooltip=["site", "wildtype"],
            color=alt.Color(
                "contact",
                scale=alt.Scale(
                    domain=["contact", "close", "distant"],
                    range=["#1f4e79", "#ff7f0e", "gray"],
                ),
                legend=alt.Legend(title="RBP Distance to Receptor"),
            ),
        )
    )
    min_ = int(df_mean[f"{metric}_E2"].min())
    max_ = int(df_mean[f"{metric}_E3"].max())
    text = (
        alt.Chart({"values": [{"x": min_, "y": max_, "text": f"r = {r_value:.2f}"}]})
        .mark_text(
            align="left",
            baseline="top",
            dx=0,  # Adjust this for position
            dy=0,  # Adjust this for position
        )
        .encode(x=alt.X("x:Q"), y=alt.Y("y:Q"), text="text:N")
    )
    plot = chart + text
    return plot


E2_E3_plot = median_correlation_plot(distance_df, "sum")
E2_E3_plot.display()
if entry_region_boxplot_plot is not None:
    E2_E3_plot.save(E2_E3_entry_corr_plot)


def correlation_plot(df):
    variant_selector = alt.selection_point(
        on="mouseover", empty=False, nearest=True, fields=["contact"], value=1
    )
    chart = (
        alt.Chart(df)
        .mark_point(opacity=0.2,size=15)
        .encode(
            x=alt.X("effect_E2", title="Cell Entry for EFNB2"),
            y=alt.Y("effect_E3", title="Cell Entry for EFNB3"),
            tooltip=["site", "wildtype", "mutant"],
            color=alt.Color(
                "contact",
                scale=alt.Scale(
                    domain=["contact", "close", "distant"],
                    range=["#1f4e79", "#ff7f0e", "gray"],
                ),
                legend=alt.Legend(title="RBP Distance to Receptor"),
            ),
            opacity=alt.condition(variant_selector,alt.value(1),alt.value(0.2)),
        ).add_params(variant_selector)
    )
    min_ = int(df["effect_E2"].min())
    max_ = int(df["effect_E3"].max())

    return chart


tmp_img_corr = correlation_plot(distance_df)
tmp_img_corr.display()
if entry_region_boxplot_plot is not None:
    tmp_img_corr.save(E2_E3_entry_all_muts_plot)

if entry_region_boxplot_plot is not None:
    (E2_E3_plot | tmp_img_corr).save(combined_E2_E3_correlation_plots)
0.81

Make boxplot showing median entry by RBP region¶

In [15]:
def make_boxplot_entry_region(df):
    barrel_ranges = {
        "Stalk": list(range(96, 147)),
        "Neck": list(range(148, 165)),
        "Linker": list(range(166, 177)),
        "Head": list(range(178, 602)),
        'Receptor Contact': config['contact_sites'],
        "Total": list(range(71, 602)),
    }
    custom_order = ["Stalk", "Neck", "Linker", "Head", "Receptor Contact", "Total"]
    empty_charts = []
    for selection in ["CHO-EFNB2", "CHO-EFNB3"]:
        if selection == "CHO-EFNB2":
            effect_name = "EFNB2"
        else:
            effect_name = "EFNB3"

        tmp_df = df[df["cell_type"] == selection]
        agg_means = []

        # For each barrel, filter the site_means dataframe to the sites belonging to that barrel and then store the means
        for barrel, sites in barrel_ranges.items():
            subset = tmp_df[tmp_df["site"].isin(sites)]
            for _, row in subset.iterrows():
                agg_means.append(
                    {"region": barrel, "effect": row["effect"], "site": row["site"]}
                )
            agg_means_df = pd.DataFrame(agg_means)

        chart = (
            alt.Chart(agg_means_df, title=f"{selection}")
            .mark_boxplot(color='darkgray')
            .encode(
                x=alt.X(
                    "region:O",
                    sort=custom_order,
                    title="RBP Region",
                    axis=alt.Axis(labelAngle=-90),
                ),
                y=alt.Y(
                    "effect",
                    title=f"Cell Entry",
                    axis=alt.Axis(grid=True, tickCount=4),
                ),
                tooltip=["region", "effect", "site"],
            ).properties(width=config['width'],height=config['height'])
        )
        empty_charts.append(chart)
    combined_effect_chart = alt.hconcat(*empty_charts).resolve_scale(
        y="shared", x="shared", color="independent"
    )
    return combined_effect_chart


entry_region_boxplot = make_boxplot_entry_region(concat_df)
entry_region_boxplot.display()
if entry_region_boxplot_plot is not None:
    entry_region_boxplot.save(entry_region_boxplot_plot)

Check for potential neutral/beneficial glycosylation sites¶

In [16]:
def find_potential_glycan_sites(df):
    filtered = df[df["wildtype"].isin(["T", "S"])]
    matching_sites = []
    for index, row in filtered.iterrows():
        # Check for existence of site two numbers prior with 'N' mutant and positive effect
        prior_rows = df[
            (df["site"] == row["site"] - 2) & (df["mutant"] == "N") & (df["effect"] > 0)
        ]
        if not prior_rows.empty:
            matching_sites.append(row["site"])
    unique_list1 = list(set(matching_sites))
    unique_list1 = [x - 2 for x in unique_list1]

    filtered = df[df["wildtype"].isin(["N"])]
    matching_sites = []
    for index, row in filtered.iterrows():
        # Check for existence of site two numbers prior with 'N' mutant and positive effect
        prior_rows = df[
            (df["site"] == row["site"] + 2)
            & (df["mutant"].isin(["T", "S"]))
            & (df["effect"] > 0)
        ]
        if not prior_rows.empty:
            matching_sites.append(row["site"])
    unique_list2 = list(set(matching_sites))
    unique_list = unique_list1 + unique_list2
    unique_list.sort()
    print(f"The total number of potential PNLG sites are: {len(unique_list)}")
    return unique_list


PNLG = find_potential_glycan_sites(func_scores_E3)

surface_df = pd.read_csv(surface)
surface_df.rename(columns={"exposure_A": "Surface Exposure"}, inplace=True)
PNLG_SURFACE = surface_df[surface_df["site"].isin(PNLG)]
PNLG_SURFACE = list(
    PNLG_SURFACE[PNLG_SURFACE["Surface Exposure"] >= 25]["site"].unique()
)

print(f"\nThe surface exposed PNLG sites are: {PNLG_SURFACE}\n")

glycans = config["glycans"]
filtered_PNLG_SURFACE = [num for num in PNLG_SURFACE if num not in glycans]

print(filtered_PNLG_SURFACE)

print(len(filtered_PNLG_SURFACE))
The total number of potential PNLG sites are: 33

The surface exposed PNLG sites are: [159, 180, 191, 192, 275, 288, 306, 311, 326, 378, 383, 403, 417, 423, 473, 478, 518, 543, 554, 570, 600]

[180, 191, 192, 275, 288, 311, 326, 383, 403, 423, 473, 478, 518, 543, 554, 570, 600]
17
In [ ]: